from scipy.stats import norm, binom
from scipy.optimize import minimize
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.mixture import GaussianMixture
# import os
# os.environ["OMP_NUM_THREADS"] = "1"作品五：探討混合常態分佈的參數估計與評估 —— MLE 與 GMM 方法之比較
作品目標：
本作品針對兩組來自不同參數的混合常態分配樣本進行實驗，旨在探討利用最大概似估計法（MLE）準確估計出兩組資料各自的參數值，並通過蒙地卡羅模擬實驗，計算參數估計值的平均值（Mean）、偏差（Bias）以及均方根誤差（RMSE），以評估 MLE 的參數估計表現。此外，本作品亦採用 sklearn.mixture.GaussianMixture 方法對資料進行參數估計，並重複相同的模擬實驗步驟，比較該方法與 MLE 在參數估計結果上的表現差異。
套件宣告：
第 1 步：混合常態參數估計
針對兩組來自不同參數的混合常態分配樣本進行實驗，探討利用最大概似估計法（MLE）準確估計出兩組資料各自的參數值。
機率密度函數（PDF）： \[ f(x|\Omega) = \pi_1 f(x | \mu_1, \sigma_1^2) + \pi_2 f(x | \mu_2, \sigma_2^2) \] 其中 \(\Omega = \{\pi_1, \pi_2, \mu_1, \sigma_1^2, \mu_2, \sigma_2^2 | \pi_1 + \pi_2 = 1\}\)
對數聯合概似函數的最大值： \[max_{\Omega = \{\pi_1, \pi_2, \mu_1, \sigma_1^2, \mu_2, \sigma_2^2 | \pi_1 + \pi_2 = 1, \pi_1, \pi_2, \sigma_1^2, \sigma_2^2 > 0 \}} L(\Omega)\] 其中目標函數為： \[ \begin{align*} L(\Omega) &= \ln \prod^{N}_{i = 1}(\pi_1 f(x_i | \mu_1, \sigma_1^2) + \pi_2 f(x_i | \mu_2, \sigma_2^2)) \\ &= \sum^{N}_{i = 1} \ln(\pi_1 f(x_i | \mu_1, \sigma_1^2) + \pi_2 f(x_i | \mu_2, \sigma_2^2)) \end{align*} \]
1.1 繪製母體參數分佈圖
分別決定出本次實驗的兩組樣本所來自的母體參數值，即混合常態分配之參數值，並繪圖呈現其分佈情況，以及加上組成該分配的兩組常態分配的分佈圖。
- 第 1 組母體參數值：\(\pi_1 = 0.8, \mu_1 = -2, \sigma_1 = 1, \mu_2 = 1, \sigma_2 = 0.5\)
- 第 2 組母體參數值：\(\pi_1 = 0.3, \mu_1 = -0.5, \sigma_1 = 1, \mu_2 = 1.5, \sigma_2 = 2\)
n = 100
pi1s = [0.8, 0.3]
mu1s, s1s = [-2, -0.5], [1, 1]
mu2s, s2s = [1, 1.5], [0.5, 2]
x_ranges = [-6, -7]
y_ranges = [4, 9]
plt.style.use('ggplot')
fig, axs = plt.subplots(1, 2, figsize = (12, 4))
for i, (pi1, mu1, s1, mu2, s2, x_range, y_range) in enumerate(zip(pi1s, mu1s, s1s, mu2s, s2s, x_ranges, y_ranges)):
f = lambda x: pi1 * norm.pdf(x, mu1, s1) + (1 - pi1) * norm.pdf(x, mu2, s2)
x = np.linspace(x_range, y_range, 1000)
y = f(x)
ax = axs.ravel()
ax[i].plot(x, y, color = '#9F88FF', label = 'Normal Mixture', linewidth = 3)
ax[i].plot(x, norm.pdf(x, mu1, s1), color = '#FF8888', linestyle = '--', linewidth = 3
, label = 'N({},{})'.format(mu1, s1))
ax[i].plot(x, norm.pdf(x, mu2, s2), color = '#99BBFF', linestyle = '--', linewidth = 3
, label = 'N({},{})'.format(mu2, s2))
ax[i].set_title('$\pi_1 = {}, \mu_1 = {}, \sigma_1 = {}, \mu_2 = {}, \sigma_2 = {}$' \
.format(pi1, mu1, s1, mu2, s2), fontsize = 12, color = 'black')
ax[i].set_xlabel('x', fontsize = 12, color = 'black')
ax[i].set_ylabel('Density', fontsize = 12, color = 'black')
ax[i].grid(True, linestyle = '-.')
ax[i].tick_params(axis = 'both', labelsize = 12, colors = 'black')
lgd = ax[i].legend(edgecolor = '#666666', prop = {'size': 11})
lgd.get_frame().set_linestyle('-.')
lgd.get_frame().set_alpha(0.4)
plt.show()
1.2 樣本數對第 1 組混合常態分佈估計的影響
生成 n 個來自混合常態的隨機樣本，其中分佈參數為 \(\pi_1 = 0.8, \mu_1 = -2, \sigma_1 = 1, \mu_2 = 1, \sigma_2 = 0.5\)，樣本數為 n = 50，100，300，500，1000，5000。針對每個樣本數 n 畫出其分佈直方圖，並利用 MLE 估計其參數值，畫出其樣本分佈圖，以及真實母體分佈圖，觀察它們之間隨著樣本數增加而產生的變化以及差異。
# 設定參數
n = [50, 100, 300, 500, 1000, 5000]
pi1 = 0.8
mu1, s1 = -2, 1
mu2, s2 = 1, 0.5
# 理論混合常態分配 pdf
f = lambda x: pi1 * norm.pdf(x, mu1, s1) + (1 - pi1) * norm.pdf(x, mu2, s2)
x = np.linspace(-6, 4, 1000)
y = f(x)
# 設定演算法初始設定
params0 = [0.8, -2, 1, 1, 0.5]
bnd = [(0, 1), (None, None), (None, None), (None, None), (None, None)]
opts = {'disp': False, 'maxiter': 8000}
plt.style.use('ggplot')
fig, axs = plt.subplots(2, 3, figsize = (12, 6))
for j in range(len(n)):
# 生成樣本
n1 = binom.rvs(n[j], pi1)
X = np.r_[norm.rvs(mu1, s1, size = n1), norm.rvs(mu2, s2, size = n[j] - n1)]
# 定義 mixed normal 對數最大概似函數
neg_log_likelihood = lambda x: -np.sum(np.log(x[0] * norm.pdf(X, x[1], x[2]) + \
(1 - x[0]) * norm.pdf(X, x[3], x[4])))
# 計算最大概似估計
result = minimize(neg_log_likelihood, params0, method = 'Nelder-Mead'
, bounds = bnd, options = opts)
ax = axs.ravel()
ax[j].plot(x, y, color = 'blue', label = 'True f(x)', linewidth = 2.5, alpha = 0.7)
ax[j].plot(x, result.x[0] * norm.pdf(x, result.x[1], result.x[2]) + \
(1 - result.x[0]) * norm.pdf(x, result.x[3], result.x[4]), label = 'Estimated f(x)'
, color = 'green', linewidth = 2.5, linestyle = '--', alpha = 0.7)
ax[j].hist(X, bins = 45, density = True, color = '#FFB347', edgecolor = 'white')
ax[j].set_title('n = {}'.format(n[j]), fontsize = 15)
ax[j].grid(True, linestyle = '-.')
ax[j].tick_params(axis = 'both', labelsize = 10, colors = 'black')
lgd = ax[j].legend(edgecolor = '#666666', prop = {'size': 10})
lgd.get_frame().set_linestyle('-.')
lgd.get_frame().set_alpha(0.4)
fig.suptitle('Population: Normal Mixture $\pi_1 = {}, \mu_1 = {}, \sigma_1 = {}, \mu_2 = {}, \sigma_2 = {}$' \
.format(pi1, mu1, s1, mu2, s2), fontsize = 17)
plt.tight_layout()
plt.show()
注意事項與討論：
從直方圖和樣本分佈的圖形可以觀察到，當樣本數 n 增加時，樣本分佈與真實分佈之間的差異逐漸縮小。
當 n = 50 時，可見直方圖中樣本分佈的雙峰特徵並不明顯，且估計曲線受隨機波動的影響較大，與真實分佈曲線之間的吻合程度不高，但其已經開始呈現出混合常態分佈的雙峰特徵。
當 n 增加至 100、300、500 時，直方圖逐漸顯示出雙峰形狀，且估計曲線與真實分佈曲線也在逐漸吻合。
當 n 達到 1000 和 5000 時，估計曲線幾乎完全貼合於真實分佈曲線，顯示出參數估計的穩定性和準確性隨著樣本數的增加而顯著提高。
最大概似估計 MLE 能夠有效估計混合常態分佈的參數，且樣本數的增加能顯著提升估計結果的準確性。
當樣本數較少時，受隨機性影響，估計的參數值可能偏離真實值，但隨著樣本數增加，參數估計會逐漸趨於穩定，並能更準確地反映真實的母體分佈。
實驗結果說明了樣本數對於參數估計的重要性，尤其在分析具有複雜結構的分佈時，大樣本的支持能更清晰地揭示分佈特徵。
1.3 樣本數對第 2 組混合常態分佈估計的影響
與前面的實驗步驟相同，但其中分佈參數改為 \(\pi_1 = 0.3, \mu_1 = -0.5, \sigma_1 = 1, \mu_2 = 1.5, \sigma_2 = 2\)。
# 設定參數
n = [50, 100, 300, 500, 1000, 5000]
pi1 = 0.3
mu1, s1 = -0.5, 1
mu2, s2 = 1.5, 2
# 理論混合常態分配 pdf
f = lambda x: pi1 * norm.pdf(x, mu1, s1) + (1 - pi1) * norm.pdf(x, mu2, s2)
x = np.linspace(-7, 9, 1000)
y = f(x)
# 設定演算法初始設定
params0 = [0.3, -0.5, 1, 1.5, 2]
bnd = [(0, 1), (None, None), (None, None), (None, None), (None, None)]
opts = {'disp': False, 'maxiter': 8000}
plt.style.use('ggplot')
fig, axs = plt.subplots(2, 3, figsize = (12, 6))
for j in range(len(n)):
# 生成樣本
n1 = binom.rvs(n[j], pi1)
X = np.r_[norm.rvs(mu1, s1, size = n1), norm.rvs(mu2, s2, size = n[j] - n1)]
# 定義 mixed normal 對數最大概似函數
neg_log_likelihood = lambda x: -np.sum(np.log(x[0] * norm.pdf(X, x[1], x[2]) + \
(1 - x[0]) * norm.pdf(X, x[3], x[4])))
# 計算最大概似估計
result = minimize(neg_log_likelihood, params0, method = 'Nelder-Mead'
, bounds = bnd, options = opts)
ax = axs.ravel()
ax[j].plot(x, y, color = '#FFA500', label = 'True f(x)', linewidth = 3)
ax[j].plot(x, result.x[0] * norm.pdf(x, result.x[1], result.x[2]) + \
(1 - result.x[0]) * norm.pdf(x, result.x[3], result.x[4]), label = 'Estimated f(x)'
, color = '#FF69B4', linestyle = '--', linewidth = 2.5)
ax[j].hist(X, bins = 45, density = True, color = '#4682B4', edgecolor = 'white')
ax[j].set_title('n = {}'.format(n[j]), fontsize = 15)
ax[j].grid(True, linestyle = '-.')
ax[j].tick_params(axis = 'both', labelsize = 10, colors = 'black')
lgd = ax[j].legend(edgecolor = '#666666', prop = {'size': 10})
lgd.get_frame().set_linestyle('-.')
lgd.get_frame().set_alpha(0.4)
fig.suptitle('Population: Normal Mixture $\pi_1 = {}, \mu_1 = {}, \sigma_1 = {}, \mu_2 = {}, \sigma_2 = {}$' \
.format(pi1, mu1, s1, mu2, s2), fontsize = 17)
plt.tight_layout()
plt.show()
注意事項與討論：
與 1.2 中的分佈相比，此次分佈的參數設置使得兩個常態分佈的重疊程度較高，分佈的雙峰特徵更加模糊。
當 n = 50 和 n = 100 時，直方圖幾乎沒有呈現出明顯的雙峰結構，難以區分兩個常態分佈的特徵。此外，估計曲線雖然已經可以看出些許的雙峰特徵，但其與真實分佈曲線的吻合程度明顯很低，容易受到隨機波動的影響。
隨著樣本數增加至 n = 300 和 n = 500，直方圖開始呈現常態分佈的基本形狀，但混合常態分佈的雙峰特徵仍然不夠明顯，且由於兩個常態分佈的重疊區域較廣，估計曲線仍略有偏差，但也開始逐漸向真實分佈曲線接近。
當樣本數進一步增加至 n = 1000 和 n = 5000 時，直方圖的形狀與真實分佈曲線逐漸吻合，估計曲線也逐漸貼合真實母體分佈，表示最大概似估計的結果幾乎完全貼合於真實分佈，但由於兩個分佈的重疊部分較大，雙峰的分界仍然不如 1.2 那麼清晰。
與 1.2 的實驗相比，當組成混合常態分佈的分佈參數之間差異較小時會使得在區分上更加困難，這也顯示了估計參數的難度會隨著分佈特性而改變。
當兩個常態分佈高度重疊時，使得小樣本下參數估計的難度更大，幾乎無法觀察到混合常態分佈的雙峰特徵，且估計參數可能會偏離真實值，進一步顯示了小樣本在高重疊分佈中的限制。
當樣本數足夠大，即 n = 5000 時，即使分佈重疊較大，MLE 仍然可以逼近真實分佈參數，但整體估計過程的穩定性和準確性明顯不如重疊較少的分佈，如 1.2。
與 1.2 的分佈相比，本次實驗更加強調了樣本數對於較為複雜的分佈結構的參數估計的重要性，尤其當混合常態分佈的組成分佈高度重疊時，大量的樣本數將會是準確描述分佈特徵的關鍵。
第 2 步：MLE 參數估計評估表現
通過蒙地卡羅模擬實驗，分別針對兩組樣本的參數估計值進行平均值（Mean）、偏差（Bias）以及均方根誤差（RMSE）的計算，評估 MLE 在參數估計上的表現。
Mean 計算：
理論值： \[ Mean(\hat{\theta}) = E(\hat{\theta}) \]
實驗值： \[ \begin{align*} Mean(\hat{\theta}) &= \frac{1}{N}\sum_{k=1}^{N}\hat{\theta}_k \\ &= \bar{\theta} \end{align*} \]
Bias 計算：
理論值： \[ Bias(\hat{\theta}) = E(\hat{\theta}) - \theta \]
實驗值： \[ \begin{align*} Bias(\hat{\theta}) &= \frac{1}{N}\sum_{k=1}^{N}\hat{\theta}_k - \theta \\ &= \bar{\theta} - \theta \end{align*} \]
Variance 計算：
理論值： \[ Var(\hat{\theta}) = E[\hat{\theta} - E(\hat{\theta})]^2 \]
實驗值： \[ Var(\hat{\theta}) = \frac{1}{N-1} \sum^{N}_{k=1} (\hat{\theta}_k - \bar{\theta})^2 \]
RMSE 計算：
理論值： \[ \begin{align*} RMSE(\hat{\theta}) &= \sqrt{MSE(\hat{\theta})} \\ &= \sqrt{E(\hat{\theta} - \theta)^2} \end{align*} \]
實驗值： \[ \begin{align*} RMSE(\hat{\theta}) &= \sqrt{MSE(\hat{\theta})} \\ &= \sqrt{\frac{1}{N} \sum^{N}_{k=1} (\hat{\theta}_k - \theta)^2} \end{align*} \]
2.1 評估不同樣本數對第 1 組分佈參數之估計表現
生成 n 個來自混合常態的隨機樣本，其中分佈參數為 \(\pi_1 = 0.8, \mu_1 = -2, \sigma_1 = 1, \mu_2 = 1, \sigma_2 = 0.5\)，樣本數為 n = 50，100，300，500，1000，5000。針對每個樣本數 n 進行 10000 次的蒙地卡羅模擬實驗，計算參數估計值的平均值（Mean）、偏差（Bias）以及均方根誤差（RMSE），以此評估 MLE 在參數估計上的表現。
# 設定參數
N = 10000 # number of simulations
n = [50, 100, 300, 500, 1000, 5000]
pi1 = 0.8
mu1, s1 = -2, 1
mu2, s2 = 1, 0.5
# 設定演算法初始設定
params0 = [0.8, -2, 1, 1, 0.5]
bnd = [(0, 1), (None, None), (None, None), (None, None), (None, None)]
opts = {'disp': False, 'maxiter': 8000}
mean_RESULT = np.zeros((len(n), len(params0)))
bias_RESULT = np.zeros((len(n), len(params0)))
rmse_RESULT = np.zeros((len(n), len(params0)))
# 添加数值稳定性处理
def safe_norm_pdf(x, loc, scale):
scale = np.clip(scale, 1e-10, None) # 防止 scale 为零或接近于零
return norm.pdf(x, loc, scale)
for j in range(len(n)):
RESULT = np.zeros((N, 5))
# 生成來自 normal mixture 的資料
n1 = binom.rvs(n[j], pi1)
SAMPLE = np.r_[norm.rvs(mu1, s1, size = (n1, N)), \
norm.rvs(mu2, s2, size = (n[j] - n1, N))]
for i in range(N):
X = SAMPLE[:, i]
neg_log_likelihood = lambda x: -np.sum(np.log(np.clip(x[0] * safe_norm_pdf(X, x[1], x[2]) + \
(1 - x[0]) * safe_norm_pdf(X, x[3], x[4]), 1e-10, None)))
result = minimize(neg_log_likelihood, params0, method = 'Nelder-Mead', \
bounds = bnd, options = opts)
RESULT[i, :] = result.x
# Evaluation
mean_RESULT[j, :] = np.mean(RESULT, axis = 0)
bias_RESULT[j, :] = mean_RESULT[j, :] - params0
rmse_RESULT[j, :] = np.sqrt(np.mean((RESULT - params0) ** 2, axis = 0))# 整理表格
columns = ['pi1={}'.format(pi1), 'mu1={}'.format(mu1), 's1={}'.format(s1)
, 'mu2={}'.format(mu2), 's2={}'.format(s2)]
columns_mean = pd.MultiIndex.from_product([['Mean for the estimated parameters at N = {}'.format(N)], \
pd.Index(columns, name = 'n')])
columns_bias = pd.MultiIndex.from_product([['Bias for the estimated parameters at N = {}'.format(N)], \
pd.Index(columns, name = 'n')])
columns_rmse = pd.MultiIndex.from_product([['RMSE for the estimated parameters at N = {}'.format(N)], \
pd.Index(columns, name = 'n')])
df_mean = pd.DataFrame(mean_RESULT, index = n, columns = columns_mean)
df_bias = pd.DataFrame(bias_RESULT, index = n, columns = columns_bias)
df_rmse = pd.DataFrame(rmse_RESULT, index = n, columns = columns_rmse)
# 使用 Styler 对象来设置列名和单元格内容的对齐方式
df_mean1 = df_mean.style.set_table_styles({
('Mean for the estimated parameters at N = {}'.format(N), 'pi1={}'.format(pi1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Mean for the estimated parameters at N = {}'.format(N), 'mu1={}'.format(mu1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Mean for the estimated parameters at N = {}'.format(N), 's1={}'.format(s1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Mean for the estimated parameters at N = {}'.format(N), 'mu2={}'.format(mu2)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Mean for the estimated parameters at N = {}'.format(N), 's2={}'.format(s2)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
}, overwrite=False).set_properties(**{'text-align': 'center'})
df_bias1 = df_bias.style.set_table_styles({
('Bias for the estimated parameters at N = {}'.format(N), 'pi1={}'.format(pi1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Bias for the estimated parameters at N = {}'.format(N), 'mu1={}'.format(mu1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Bias for the estimated parameters at N = {}'.format(N), 's1={}'.format(s1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Bias for the estimated parameters at N = {}'.format(N), 'mu2={}'.format(mu2)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Bias for the estimated parameters at N = {}'.format(N), 's2={}'.format(s2)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
}, overwrite=False).set_properties(**{'text-align': 'center'})
df_rmse1 = df_rmse.style.set_table_styles({
('RMSE for the estimated parameters at N = {}'.format(N), 'pi1={}'.format(pi1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('RMSE for the estimated parameters at N = {}'.format(N), 'mu1={}'.format(mu1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('RMSE for the estimated parameters at N = {}'.format(N), 's1={}'.format(s1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('RMSE for the estimated parameters at N = {}'.format(N), 'mu2={}'.format(mu2)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('RMSE for the estimated parameters at N = {}'.format(N), 's2={}'.format(s2)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
}, overwrite=False).set_properties(**{'text-align': 'center'})
# 添加水平線在每個 row
df_mean1.set_table_styles([{'selector': 'tr', 'props': [('border-top', '1px solid #666666')]}], overwrite=False)
df_bias1.set_table_styles([{'selector': 'tr', 'props': [('border-top', '1px solid #666666')]}], overwrite=False)
df_rmse1.set_table_styles([{'selector': 'tr', 'props': [('border-top', '1px solid #666666')]}], overwrite=False)
df_mean1.set_table_styles([{'selector': 'tr', 'props': [('border-bottom', '1px solid #666666')]}], overwrite=False)
df_bias1.set_table_styles([{'selector': 'tr', 'props': [('border-bottom', '1px solid #666666')]}], overwrite=False)
df_rmse1.set_table_styles([{'selector': 'tr', 'props': [('border-bottom', '1px solid #666666')]}], overwrite=False)
# 显示 DataFrame
display(df_mean1)
display(df_bias1)
display(df_rmse1)| Mean for the estimated parameters at N = 10000 | |||||
|---|---|---|---|---|---|
| n | pi1=0.8 | mu1=-2 | s1=1 | mu2=1 | s2=0.5 |
| 50 | 0.773740 | -2.009000 | 0.967260 | 0.979545 | 0.465184 |
| 100 | 0.826820 | -2.006365 | 0.986970 | 0.985783 | 0.478259 |
| 300 | 0.779784 | -1.999313 | 0.996391 | 0.998326 | 0.493434 |
| 500 | 0.799751 | -1.999851 | 0.998492 | 0.998532 | 0.496584 |
| 1000 | 0.799870 | -2.000500 | 0.998825 | 0.998953 | 0.498085 |
| 5000 | 0.798599 | -1.999918 | 0.999699 | 0.999793 | 0.499637 |
| Bias for the estimated parameters at N = 10000 | |||||
|---|---|---|---|---|---|
| n | pi1=0.8 | mu1=-2 | s1=1 | mu2=1 | s2=0.5 |
| 50 | -0.026260 | -0.009000 | -0.032740 | -0.020455 | -0.034816 |
| 100 | 0.026820 | -0.006365 | -0.013030 | -0.014217 | -0.021741 |
| 300 | -0.020216 | 0.000687 | -0.003609 | -0.001674 | -0.006566 |
| 500 | -0.000249 | 0.000149 | -0.001508 | -0.001468 | -0.003416 |
| 1000 | -0.000130 | -0.000500 | -0.001175 | -0.001047 | -0.001915 |
| 5000 | -0.001401 | 0.000082 | -0.000301 | -0.000207 | -0.000363 |
| RMSE for the estimated parameters at N = 10000 | |||||
|---|---|---|---|---|---|
| n | pi1=0.8 | mu1=-2 | s1=1 | mu2=1 | s2=0.5 |
| 50 | 0.061016 | 0.205564 | 0.160085 | 0.252441 | 0.175287 |
| 100 | 0.043892 | 0.137564 | 0.105532 | 0.199863 | 0.140579 |
| 300 | 0.025060 | 0.077644 | 0.061129 | 0.083469 | 0.061001 |
| 500 | 0.011113 | 0.059082 | 0.046000 | 0.066979 | 0.049796 |
| 1000 | 0.007671 | 0.041455 | 0.032520 | 0.047317 | 0.035387 |
| 5000 | 0.003708 | 0.018522 | 0.014677 | 0.021012 | 0.015571 |
注意事項與討論：
從表中可見，混合常態分佈參數的最大概似估計 MLE 在不同樣本數 n 下，其平均值皆相差不大，無論樣本數的大小，估計出來的參數值都很接近真實值。
但是，從偏差和 RMSE 的表中則可以明顯看出隨著樣本數的增加，其值皆會隨之下降，表示參數估計會受到樣本數的影響。
本次實驗表明，混合常態分佈的 MLE 參數估計確實會受樣本數影響。在小樣本下，參數估計結果較不穩定，偏差和 RMSE 較高；而在樣本數增加後，估計的穩定性和準確性均有顯著提升，進一步驗證了前一目標的結論。
此外，本次實驗顯示，無論樣本數大小，各參數的平均值皆接近母體真實值。這是因為組成混合常態的兩個分佈參數差異較大，即使在小樣本下，只要模擬次數足夠，其參數估計仍能準確區分這兩個分佈，因此平均值在小樣本下也得以接近真實值。
整體而言，MLE 在混合常態分佈的參數估計中表現穩定且準確，而樣本數的增加則是降低偏差與 RMSE 的關鍵。因此，進行參數估計時，建議採用較大的樣本數以提升準確性。
2.2 評估不同樣本數對第 2 組分佈參數之估計表現
與前面的實驗步驟相同，但其中分佈參數改為 \(\pi_1 = 0.3, \mu_1 = -0.5, \sigma_1 = 1, \mu_2 = 1.5, \sigma_2 = 2\)。
# 設定參數
N = 10000 # number of simulations
n = [50, 100, 300, 500, 1000, 5000]
pi1 = 0.3
mu1, s1 = -0.5, 1
mu2, s2 = 1.5, 2
# 設定演算法初始設定
params0 = [0.3, -0.5, 1, 1.5, 2]
bnd = [(0, 1), (None, None), (None, None), (None, None), (None, None)]
opts = {'disp': False, 'maxiter': 8000}
mean_RESULT = np.zeros((len(n), len(params0)))
bias_RESULT = np.zeros((len(n), len(params0)))
rmse_RESULT = np.zeros((len(n), len(params0)))
# 添加数值稳定性处理
def safe_norm_pdf(x, loc, scale):
scale = np.clip(scale, 1e-10, None) # 防止 scale 为零或接近于零
return norm.pdf(x, loc, scale)
for j in range(len(n)):
RESULT = np.zeros((N, 5))
# 生成來自 normal mixture 的資料
n1 = binom.rvs(n[j], pi1)
SAMPLE = np.r_[norm.rvs(mu1, s1, size = (n1, N)), \
norm.rvs(mu2, s2, size = (n[j] - n1, N))]
for i in range(N):
X = SAMPLE[:, i]
# neg_log_likelihood = lambda x: -np.sum(np.log(x[0] * norm.pdf(X, x[1], x[2]) + \
# (1 - x[0]) * norm.pdf(X, x[3], x[4])))
neg_log_likelihood = lambda x: -np.sum(np.log(np.clip(x[0] * safe_norm_pdf(X, x[1], x[2]) + \
(1 - x[0]) * safe_norm_pdf(X, x[3], x[4]), 1e-10, None)))
result = minimize(neg_log_likelihood, params0, method = 'Nelder-Mead', \
bounds = bnd, options = opts)
RESULT[i, :] = result.x
# Evaluation
mean_RESULT[j, :] = np.mean(RESULT, axis = 0)
bias_RESULT[j, :] = mean_RESULT[j, :] - params0
rmse_RESULT[j, :] = np.sqrt(np.mean((RESULT - params0) ** 2, axis = 0))# 整理表格
columns = ['pi1={}'.format(pi1), 'mu1={}'.format(mu1), 's1={}'.format(s1)
, 'mu2={}'.format(mu2), 's2={}'.format(s2)]
columns_mean = pd.MultiIndex.from_product([['Mean for the estimated parameters at N = {}'.format(N)], \
pd.Index(columns, name = 'n')])
columns_bias = pd.MultiIndex.from_product([['Bias for the estimated parameters at N = {}'.format(N)], \
pd.Index(columns, name = 'n')])
columns_rmse = pd.MultiIndex.from_product([['RMSE for the estimated parameters at N = {}'.format(N)], \
pd.Index(columns, name = 'n')])
df_mean = pd.DataFrame(mean_RESULT, index = n, columns = columns_mean)
df_bias = pd.DataFrame(bias_RESULT, index = n, columns = columns_bias)
df_rmse = pd.DataFrame(rmse_RESULT, index = n, columns = columns_rmse)
# 使用 Styler 对象来设置列名和单元格内容的对齐方式
df_mean1 = df_mean.style.set_table_styles({
('Mean for the estimated parameters at N = {}'.format(N), 'pi1={}'.format(pi1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Mean for the estimated parameters at N = {}'.format(N), 'mu1={}'.format(mu1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Mean for the estimated parameters at N = {}'.format(N), 's1={}'.format(s1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Mean for the estimated parameters at N = {}'.format(N), 'mu2={}'.format(mu2)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Mean for the estimated parameters at N = {}'.format(N), 's2={}'.format(s2)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
}, overwrite=False).set_properties(**{'text-align': 'center'})
df_bias1 = df_bias.style.set_table_styles({
('Bias for the estimated parameters at N = {}'.format(N), 'pi1={}'.format(pi1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Bias for the estimated parameters at N = {}'.format(N), 'mu1={}'.format(mu1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Bias for the estimated parameters at N = {}'.format(N), 's1={}'.format(s1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Bias for the estimated parameters at N = {}'.format(N), 'mu2={}'.format(mu2)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Bias for the estimated parameters at N = {}'.format(N), 's2={}'.format(s2)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
}, overwrite=False).set_properties(**{'text-align': 'center'})
df_rmse1 = df_rmse.style.set_table_styles({
('RMSE for the estimated parameters at N = {}'.format(N), 'pi1={}'.format(pi1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('RMSE for the estimated parameters at N = {}'.format(N), 'mu1={}'.format(mu1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('RMSE for the estimated parameters at N = {}'.format(N), 's1={}'.format(s1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('RMSE for the estimated parameters at N = {}'.format(N), 'mu2={}'.format(mu2)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('RMSE for the estimated parameters at N = {}'.format(N), 's2={}'.format(s2)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
}, overwrite=False).set_properties(**{'text-align': 'center'})
# 添加水平線在每個 row
df_mean1.set_table_styles([{'selector': 'tr', 'props': [('border-top', '1px solid #666666')]}], overwrite=False)
df_bias1.set_table_styles([{'selector': 'tr', 'props': [('border-top', '1px solid #666666')]}], overwrite=False)
df_rmse1.set_table_styles([{'selector': 'tr', 'props': [('border-top', '1px solid #666666')]}], overwrite=False)
df_mean1.set_table_styles([{'selector': 'tr', 'props': [('border-bottom', '1px solid #666666')]}], overwrite=False)
df_bias1.set_table_styles([{'selector': 'tr', 'props': [('border-bottom', '1px solid #666666')]}], overwrite=False)
df_rmse1.set_table_styles([{'selector': 'tr', 'props': [('border-bottom', '1px solid #666666')]}], overwrite=False)
# 显示 DataFrame
display(df_mean1)
display(df_bias1)
display(df_rmse1)| Mean for the estimated parameters at N = 10000 | |||||
|---|---|---|---|---|---|
| n | pi1=0.3 | mu1=-0.5 | s1=1 | mu2=1.5 | s2=2 |
| 50 | 0.457214 | -0.392147 | 0.964834 | 2.524397 | 1.499656 |
| 100 | 0.476407 | -0.367772 | 0.973800 | 2.238184 | 1.680808 |
| 300 | 0.396133 | -0.409818 | 0.999474 | 1.884593 | 1.856753 |
| 500 | 0.359095 | -0.437408 | 1.006318 | 1.754697 | 1.907502 |
| 1000 | 0.331196 | -0.463092 | 1.012992 | 1.645061 | 1.949769 |
| 5000 | 0.294702 | -0.496174 | 1.002232 | 1.518742 | 1.994619 |
| Bias for the estimated parameters at N = 10000 | |||||
|---|---|---|---|---|---|
| n | pi1=0.3 | mu1=-0.5 | s1=1 | mu2=1.5 | s2=2 |
| 50 | 0.157214 | 0.107853 | -0.035166 | 1.024397 | -0.500344 |
| 100 | 0.176407 | 0.132228 | -0.026200 | 0.738184 | -0.319192 |
| 300 | 0.096133 | 0.090182 | -0.000526 | 0.384593 | -0.143247 |
| 500 | 0.059095 | 0.062592 | 0.006318 | 0.254697 | -0.092498 |
| 1000 | 0.031196 | 0.036908 | 0.012992 | 0.145061 | -0.050231 |
| 5000 | -0.005298 | 0.003826 | 0.002232 | 0.018742 | -0.005381 |
| RMSE for the estimated parameters at N = 10000 | |||||
|---|---|---|---|---|---|
| n | pi1=0.3 | mu1=-0.5 | s1=1 | mu2=1.5 | s2=2 |
| 50 | 0.325228 | 0.891936 | 0.518010 | 1.731392 | 0.798179 |
| 100 | 0.302563 | 0.540694 | 0.391532 | 1.410807 | 0.599995 |
| 300 | 0.214678 | 0.341806 | 0.286134 | 0.853977 | 0.326285 |
| 500 | 0.170352 | 0.265562 | 0.234744 | 0.633321 | 0.231810 |
| 1000 | 0.120812 | 0.182444 | 0.167767 | 0.421792 | 0.151667 |
| 5000 | 0.039021 | 0.068541 | 0.069282 | 0.112181 | 0.039595 |
注意事項與討論：
從表中可見，與 2.1 中的結果相比，本次實驗的參數估計整體表現較差，此次可見在小樣本情況下，偏差和 RMSE 明顯高於 2.1。
當樣本數為 n = 50 和 n = 100 時，參數的偏差明顯較大，平均值明顯偏離母體真實值，反映出估計的穩定性較差。這是因為兩個常態分佈的重疊區域較廣，導致混合常態分佈的特徵難以被準確捕捉。
隨著樣本數增加至 n = 300、500 和 1000，估計結果的穩定性逐步提升，偏差和 RMSE 明顯逐漸下降，但整體平均值仍與真實值存在一定差異，表明高重疊分佈下的估計難度仍然存在。
當樣本數達到 n = 5000 時，參數估計的平均值接近母體真實值，偏差和 RMSE 明顯下降，但由於分佈的雙峰結構仍然不夠明顯，其偏差和 RMSE 仍不如 2.1 來得低。
與 2.1 的結果相比，本次實驗中組成混合常態分佈的兩個組成分佈參數差異較小，導致重疊程度較高，為 MLE 的參數估計增加了難度。在小樣本情況下，估計值的偏差和 RMSE 顯著增加，反映出高度重疊分佈下參數估計的穩定性較低。
隨著樣本數的增加，MLE 的估計性能逐步提升，但整體的估計準確性仍受分佈特徵影響較大。僅在樣本數達到 n = 5000 時，估計值的平均值才接近母體真實值，顯示出高度重疊分佈需要更大的樣本數以提升估計準確性。
本次實驗進一步證實，當混合常態分佈的組成分佈參數差異較小時，樣本數對於降低偏差和 RMSE 的重要性尤為重要。因此，在使用 MLE 進行參數估計時，需特別注意分佈參數的特徵，並確保樣本數充足，以提升模型的估計穩定性與準確性。
第 3 步：GMM 參數估計評估表現
採用 sklearn.mixture.GaussianMixture（後續簡稱 GMM）方法對資料進行參數估計，並重複前面兩個目標的模擬實驗步驟，比較該方法與 MLE 在參數估計結果上的表現差異。
3.1 比較 MLE 與 GMM 的參數估計
與第 1 步的 1.2 實驗步驟相同，生成 n 個來自混合常態的隨機樣本，其中分佈參數為 \(\pi_1 = 0.8, \mu_1 = -2, \sigma_1 = 1, \mu_2 = 1, \sigma_2 = 0.5\)，樣本數為 n = 50，100，300，500，1000，5000。針對每個樣本數 n 畫出其分佈直方圖，並利用 MLE 估計其參數值，畫出其樣本分佈圖，以及真實母體分佈圖。此外，利用 GMM 估計出參數值，並一樣畫出其分佈圖，觀察 MLE 與 GMM 在參數估計上的差異。
# 設定參數
n = [50, 100, 300, 500, 1000, 5000]
pi1 = 0.8
mu1, s1 = -2, 1
mu2, s2 = 1, 0.5
# 理論混合常態分配 pdf
f = lambda x: pi1 * norm.pdf(x, mu1, s1) + (1 - pi1) * norm.pdf(x, mu2, s2)
x = np.linspace(-6, 4, 1000)
y = f(x)
# 設定演算法初始設定
params0 = [0.8, -2, 1, 1, 0.5]
bnd = [(0, 1), (None, None), (None, None), (None, None), (None, None)]
opts = {'disp': False, 'maxiter': 8000}
plt.style.use('ggplot')
fig, axs = plt.subplots(2, 3, figsize = (14, 6))
for j in range(len(n)):
# 生成樣本
n1 = binom.rvs(n[j], pi1)
X = np.r_[norm.rvs(mu1, s1, size = n1), norm.rvs(mu2, s2, size = n[j] - n1)]
# 定義 mixed normal 對數最大概似函數
neg_log_likelihood = lambda x: -np.sum(np.log(x[0] * norm.pdf(X, x[1], x[2]) + \
(1 - x[0]) * norm.pdf(X, x[3], x[4])))
# 計算最大概似估計
result = minimize(neg_log_likelihood, params0, method = 'Nelder-Mead'
, bounds = bnd, options = opts)
# Gaussian Mixture 演算法
gmm = GaussianMixture(n_components = 2, covariance_type = 'full', \
verbose = 0, max_iter = 1000, tol = 1e-6)
gmm.fit(X.reshape(-1, 1)) # sample is the normal mixture data as a column vector
ax = axs.ravel()
ax[j].plot(x, y, color = '#6c7b8b', label = 'True f(x)', linewidth = 3)
ax[j].plot(x, result.x[0] * norm.pdf(x, result.x[1], result.x[2]) + \
(1 - result.x[0]) * norm.pdf(x, result.x[3], result.x[4]), label = 'Optimise Minimise'
, color = '#f0e68c', linestyle = '--', linewidth = 2.5)
ax[j].plot(x, gmm.weights_[0] * norm.pdf(x, gmm.means_[0][0], np.sqrt(gmm.covariances_[0][0])) + \
gmm.weights_[1] * norm.pdf(x, gmm.means_[1][0], np.sqrt(gmm.covariances_[1][0])) \
, label = 'Gaussian Mixture', color = '#9caf88', linestyle = '-.', linewidth = 2.5)
ax[j].hist(X, bins = 45, density = True, color = '#b0c4de', edgecolor = 'white')
ax[j].set_title('n = {}'.format(n[j]), fontsize = 15)
ax[j].grid(True, linestyle = '-.')
ax[j].tick_params(axis = 'both', labelsize = 10, colors = 'black')
lgd = ax[j].legend(edgecolor = '#666666', prop = {'size': 10})
lgd.get_frame().set_linestyle('-.')
lgd.get_frame().set_alpha(0.4)
fig.suptitle('Population: Normal Mixture $\pi_1 = {}, \mu_1 = {}, \sigma_1 = {}, \mu_2 = {}, \sigma_2 = {}$' \
.format(pi1, mu1, s1, mu2, s2), fontsize = 17)
plt.tight_layout()
plt.show()
注意事項與討論：
從圖中可以看出，無論樣本數的大小，MLE 和 GMM 的估計曲線均顯示出高度重疊的情況。
此外，兩者在樣本數較少時，估計曲線與真實分佈曲線之間存在明顯差異，但隨著樣本數的增加，估計曲線逐漸趨近真實分佈曲線。
根據實驗結果，未能顯著區分 MLE 和 GMM 之間的差異，這可能是因為混合常態分佈的組成參數差異較大，從而使得兩者在區分分佈時呈現相似的結果。
由此可見，GMM 在參數估計方面同樣表現出色，尤其是在分佈參數差異較大且樣本數較多的情況下，兩種方法均能準確地還原母體分佈的特徵，包括雙峰結構和分佈形狀。
3.2 評估 MLE 與 GMM 的參數估計表現
生成 n 個來自混合常態的隨機樣本，其中分佈參數為 \(\pi_1 = 0.8, \mu_1 = -2, \sigma_1 = 1, \mu_2 = 1, \sigma_2 = 0.5\)，樣本數為 n = 50，100，300，500，1000，5000。針對每個樣本數 n 進行 10000 次的蒙地卡羅模擬實驗，並分別使用 MLE 和 GMM 去估計參數，計算各自參數估計值的平均值（Mean）、偏差（Bias）以及均方根誤差（RMSE），以此評估 MLE 和 GMM 在參數估計上的表現。
# 設定參數
N = 10000 # number of simulations
n = [50, 100, 300, 500, 1000, 5000]
pi1 = 0.8
mu1, s1 = -2, 1
mu2, s2 = 1, 0.5
# 設定演算法初始設定
params0 = [0.8, -2, 1, 1, 0.5]
bnd = [(0, 1), (None, None), (None, None), (None, None), (None, None)]
opts = {'disp': False, 'maxiter': 8000}
mean_mini = np.zeros((len(n), len(params0))); mean_gmm = np.zeros((len(n), len(params0)))
bias_mini = np.zeros((len(n), len(params0))); bias_gmm = np.zeros((len(n), len(params0)))
rmse_mini = np.zeros((len(n), len(params0))); rmse_gmm = np.zeros((len(n), len(params0)))
for j in range(len(n)):
# 生成結果矩陣
RESULT_mini = np.zeros((N, 5))
RESULT_gmm = np.zeros((N, 5))
# 生成來自 mix beta distribution 的資料
n1 = binom.rvs(n[j], pi1)
SAMPLE = np.r_[norm.rvs(mu1, s1, size = (n1, N)), \
norm.rvs(mu2, s2, size = (n[j] - n1, N))]
for i in range(N):
X = SAMPLE[:, i]
neg_log_likelihood = lambda x: -np.sum(np.log(x[0] * norm.pdf(X, x[1], x[2]) + \
(1 - x[0]) * norm.pdf(X, x[3], x[4])))
result = minimize(neg_log_likelihood, params0, method = 'Nelder-Mead', \
bounds = bnd, options = opts)
RESULT_mini[i, :] = result.x
# Gaussian Mixture 演算法
gmm = GaussianMixture(n_components = 2, covariance_type = 'full', \
weights_init = [pi1, 1-pi1], \
means_init = np.array([[mu1], [mu2]]).reshape(2, 1), \
verbose = 0, max_iter = 8000, tol = 1e-6)
gmm.fit(X.reshape(-1, 1))
a, b = np.argmax(gmm.weights_), np.argmin(gmm.weights_)
RESULT_gmm[i, :] = [gmm.weights_[a], gmm.means_[0][0], np.sqrt(gmm.covariances_[0][0][0]), \
gmm.means_[1][0], np.sqrt(gmm.covariances_[1][0][0])]
# minimize
mean_mini[j, :] = np.mean(RESULT_mini, axis = 0)
bias_mini[j, :] = mean_mini[j, :] - params0
rmse_mini[j, :] = np.sqrt(np.mean((RESULT_mini - params0) ** 2, axis = 0))
# gmm
mean_gmm[j, :] = np.mean(RESULT_gmm, axis = 0)
bias_gmm[j, :] = mean_gmm[j, :] - params0
rmse_gmm[j, :] = np.sqrt(np.mean((RESULT_gmm - params0) ** 2, axis = 0))# 整理表格
columns = ['pi1={}'.format(pi1), 'mu1={}'.format(mu1), 's1={}'.format(s1)
, 'mu2={}'.format(mu2), 's2={}'.format(s2)]
# 创建 MultiIndex 列名
columns_mean = pd.MultiIndex.from_product([['Mean for the estimated parameters at N = {}'.format(N)], pd.Index(columns, name='n')])
columns_bias = pd.MultiIndex.from_product([['Bias for the estimated parameters at N = {}'.format(N)], pd.Index(columns, name='n')])
columns_rmse = pd.MultiIndex.from_product([['RMSE for the estimated parameters at N = {}'.format(N)], pd.Index(columns, name='n')])
# 创建新的 MultiIndex
methods = np.repeat(['optimize.minimize', 'GaussianMixture'], len(n))
ns = np.tile(n, 2)
multi_index = pd.MultiIndex.from_arrays([methods, ns])
# 创建 DataFrame 并设置 MultiIndex
df_mean = pd.DataFrame(np.vstack([mean_mini, mean_gmm]), index=multi_index, columns=columns_mean)
df_bias = pd.DataFrame(np.vstack([bias_mini, bias_gmm]), index=multi_index, columns=columns_bias)
df_rmse = pd.DataFrame(np.vstack([rmse_mini, rmse_gmm]), index=multi_index, columns=columns_rmse)
# 使用 Styler 对象来设置列名和单元格内容的对齐方式
df_mean1 = df_mean.style.set_table_styles({
('Mean for the estimated parameters at N = {}'.format(N), 'pi1={}'.format(pi1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Mean for the estimated parameters at N = {}'.format(N), 'mu1={}'.format(mu1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Mean for the estimated parameters at N = {}'.format(N), 's1={}'.format(s1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Mean for the estimated parameters at N = {}'.format(N), 'mu2={}'.format(mu2)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Mean for the estimated parameters at N = {}'.format(N), 's2={}'.format(s2)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
}, overwrite=False).set_properties(**{'text-align': 'center'})
df_bias1 = df_bias.style.set_table_styles({
('Bias for the estimated parameters at N = {}'.format(N), 'pi1={}'.format(pi1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Bias for the estimated parameters at N = {}'.format(N), 'mu1={}'.format(mu1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Bias for the estimated parameters at N = {}'.format(N), 's1={}'.format(s1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Bias for the estimated parameters at N = {}'.format(N), 'mu2={}'.format(mu2)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Bias for the estimated parameters at N = {}'.format(N), 's2={}'.format(s2)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
}, overwrite=False).set_properties(**{'text-align': 'center'})
df_rmse1 = df_rmse.style.set_table_styles({
('RMSE for the estimated parameters at N = {}'.format(N), 'pi1={}'.format(pi1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('RMSE for the estimated parameters at N = {}'.format(N), 'mu1={}'.format(mu1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('RMSE for the estimated parameters at N = {}'.format(N), 's1={}'.format(s1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('RMSE for the estimated parameters at N = {}'.format(N), 'mu2={}'.format(mu2)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('RMSE for the estimated parameters at N = {}'.format(N), 's2={}'.format(s2)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
}, overwrite=False).set_properties(**{'text-align': 'center'})
# 添加水平線在每個 row
df_mean1.set_table_styles([{'selector': 'tr', 'props': [('border-top', '1px solid #666666')]}], overwrite=False)
df_bias1.set_table_styles([{'selector': 'tr', 'props': [('border-top', '1px solid #666666')]}], overwrite=False)
df_rmse1.set_table_styles([{'selector': 'tr', 'props': [('border-top', '1px solid #666666')]}], overwrite=False)
df_mean1.set_table_styles([{'selector': 'tr', 'props': [('border-bottom', '1px solid #666666')]}], overwrite=False)
df_bias1.set_table_styles([{'selector': 'tr', 'props': [('border-bottom', '1px solid #666666')]}], overwrite=False)
df_rmse1.set_table_styles([{'selector': 'tr', 'props': [('border-bottom', '1px solid #666666')]}], overwrite=False)
# 显示 DataFrame
display(df_mean1)
display(df_bias1)
display(df_rmse1)| Mean for the estimated parameters at N = 10000 | ||||||
|---|---|---|---|---|---|---|
| n | pi1=0.8 | mu1=-2 | s1=1 | mu2=1 | s2=0.5 | |
| optimize.minimize | 50 | 0.754867 | -2.010923 | 0.969151 | 0.980275 | 0.469128 |
| 100 | 0.825918 | -2.007108 | 0.983263 | 0.986594 | 0.480410 | |
| 300 | 0.766490 | -1.999965 | 0.996626 | 0.998892 | 0.494705 | |
| 500 | 0.791774 | -2.000591 | 0.997823 | 0.998199 | 0.496855 | |
| 1000 | 0.780913 | -1.999648 | 0.999331 | 1.000191 | 0.498479 | |
| 5000 | 0.793792 | -2.000004 | 0.999834 | 1.000161 | 0.499760 | |
| GaussianMixture | 50 | 0.754222 | -2.014731 | 0.965907 | 0.973527 | 0.474352 |
| 100 | 0.825378 | -2.009128 | 0.981726 | 0.982695 | 0.483341 | |
| 300 | 0.766207 | -2.000764 | 0.995920 | 0.997906 | 0.495441 | |
| 500 | 0.791487 | -2.001375 | 0.997122 | 0.997079 | 0.497681 | |
| 1000 | 0.780640 | -2.000408 | 0.998647 | 0.999180 | 0.499232 | |
| 5000 | 0.793527 | -2.000734 | 0.999170 | 0.999120 | 0.500531 | |
| Bias for the estimated parameters at N = 10000 | ||||||
|---|---|---|---|---|---|---|
| n | pi1=0.8 | mu1=-2 | s1=1 | mu2=1 | s2=0.5 | |
| optimize.minimize | 50 | -0.045133 | -0.010923 | -0.030849 | -0.019725 | -0.030872 |
| 100 | 0.025918 | -0.007108 | -0.016737 | -0.013406 | -0.019590 | |
| 300 | -0.033510 | 0.000035 | -0.003374 | -0.001108 | -0.005295 | |
| 500 | -0.008226 | -0.000591 | -0.002177 | -0.001801 | -0.003145 | |
| 1000 | -0.019087 | 0.000352 | -0.000669 | 0.000191 | -0.001521 | |
| 5000 | -0.006208 | -0.000004 | -0.000166 | 0.000161 | -0.000240 | |
| GaussianMixture | 50 | -0.045778 | -0.014731 | -0.034093 | -0.026473 | -0.025648 |
| 100 | 0.025378 | -0.009128 | -0.018274 | -0.017305 | -0.016659 | |
| 300 | -0.033793 | -0.000764 | -0.004080 | -0.002094 | -0.004559 | |
| 500 | -0.008513 | -0.001375 | -0.002878 | -0.002921 | -0.002319 | |
| 1000 | -0.019360 | -0.000408 | -0.001353 | -0.000820 | -0.000768 | |
| 5000 | -0.006473 | -0.000734 | -0.000830 | -0.000880 | 0.000531 | |
| RMSE for the estimated parameters at N = 10000 | ||||||
|---|---|---|---|---|---|---|
| n | pi1=0.8 | mu1=-2 | s1=1 | mu2=1 | s2=0.5 | |
| optimize.minimize | 50 | 0.068595 | 0.206615 | 0.159463 | 0.230084 | 0.164913 |
| 100 | 0.045785 | 0.137913 | 0.109525 | 0.208402 | 0.142789 | |
| 300 | 0.036676 | 0.078287 | 0.061762 | 0.080085 | 0.059925 | |
| 500 | 0.013818 | 0.059226 | 0.045940 | 0.065333 | 0.049125 | |
| 1000 | 0.020614 | 0.041754 | 0.033532 | 0.044741 | 0.033135 | |
| 5000 | 0.007073 | 0.018562 | 0.014696 | 0.020314 | 0.015275 | |
| GaussianMixture | 50 | 0.066018 | 0.205685 | 0.159985 | 0.238590 | 0.169569 |
| 100 | 0.043398 | 0.137470 | 0.108964 | 0.209344 | 0.143391 | |
| 300 | 0.036924 | 0.078117 | 0.061609 | 0.080200 | 0.060023 | |
| 500 | 0.014004 | 0.059164 | 0.045875 | 0.065580 | 0.049314 | |
| 1000 | 0.020870 | 0.041695 | 0.033473 | 0.044879 | 0.033263 | |
| 5000 | 0.007309 | 0.018550 | 0.014686 | 0.020393 | 0.015357 | |
注意事項與討論：
從結果中可以看到，與 3.1 的情況相似，無論樣本數的大小，MLE 和 GMM 在參數估計上的表現幾乎沒有顯著差異。
此外，兩者在樣本數較少時，估計的平均值與真實值之間仍然存在差異，但隨著樣本數的增加，平均值會逐漸接近真實值，偏差和 RMSE 也隨之降低。
實驗結果與 3.1 類似，未能明顯顯示 MLE 和 GMM 之間的差異。原因與 3.1 相同，兩個常態分佈的參數差異較大，因此能夠較容易區分，無法顯現兩者在處理分佈細節上的不同。
兩者在樣本數增大後，都展現出估計準確性的提高，並且在大樣本下，其估計結果與真實母體參數值相當接近。因此，GMM 同樣是一個可靠的參數估計方法。
3.3 比較 MLE 與 GMM 的參數估計
與 3.1 的實驗步驟相同，但其中分佈參數改為 \(\pi_1 = 0.3, \mu_1 = -0.5, \sigma_1 = 1, \mu_2 = 1.5, \sigma_2 = 2\)，觀察 MLE 與 GMM 在參數估計上的差異。
# 設定參數
n = [50, 100, 300, 500, 1000, 5000]
pi1 = 0.3
mu1, s1 = -0.5, 1
mu2, s2 = 1.5, 2
# 理論混合常態分配 pdf
f = lambda x: pi1 * norm.pdf(x, mu1, s1) + (1 - pi1) * norm.pdf(x, mu2, s2)
x = np.linspace(-7, 9, 1000)
y = f(x)
# 設定演算法初始設定
params0 = [0.3, -0.5, 1, 1.5, 2]
bnd = [(0, 1), (None, None), (None, None), (None, None), (None, None)]
opts = {'disp': False, 'maxiter': 8000}
plt.style.use('ggplot')
fig, axs = plt.subplots(2, 3, figsize = (14, 6))
for j in range(len(n)):
# 生成樣本
n1 = binom.rvs(n[j], pi1)
X = np.r_[norm.rvs(mu1, s1, size = n1), norm.rvs(mu2, s2, size = n[j] - n1)]
# 定義 mixed normal 對數最大概似函數
neg_log_likelihood = lambda x: -np.sum(np.log(x[0] * norm.pdf(X, x[1], x[2]) + \
(1 - x[0]) * norm.pdf(X, x[3], x[4])))
# 計算最大概似估計
result = minimize(neg_log_likelihood, params0, method = 'Nelder-Mead'
, bounds = bnd, options = opts)
# Gaussian Mixture 演算法
gmm = GaussianMixture(n_components = 2, covariance_type = 'full', \
verbose = 0, max_iter = 1000, tol = 1e-6)
gmm.fit(X.reshape(-1, 1)) # sample is the normal mixture data as a column vector
ax = axs.ravel()
ax[j].plot(x, y, color = '#d62728', label = 'True f(x)', linewidth = 3)
ax[j].plot(x, result.x[0] * norm.pdf(x, result.x[1], result.x[2]) + \
(1 - result.x[0]) * norm.pdf(x, result.x[3], result.x[4]), label = 'Optimise Minimise'
, color = '#ffbf00', linestyle = '--', linewidth = 2.5)
ax[j].plot(x, gmm.weights_[0] * norm.pdf(x, gmm.means_[0][0], np.sqrt(gmm.covariances_[0][0])) + \
gmm.weights_[1] * norm.pdf(x, gmm.means_[1][0], np.sqrt(gmm.covariances_[1][0])) \
, label = 'Gaussian Mixture', color = '#17becf', linestyle = '-.', linewidth = 2.5)
ax[j].hist(X, bins = 45, density = True, color = '#4d4d4d', edgecolor = 'white', alpha = 0.5)
ax[j].set_title('n = {}'.format(n[j]), fontsize = 15)
ax[j].grid(True, linestyle = '-.')
ax[j].tick_params(axis = 'both', labelsize = 10, colors = 'black')
lgd = ax[j].legend(edgecolor = '#666666', prop = {'size': 10})
lgd.get_frame().set_linestyle('-.')
lgd.get_frame().set_alpha(0.4)
fig.suptitle('Population: Normal Mixture $\pi_1 = {}, \mu_1 = {}, \sigma_1 = {}, \mu_2 = {}, \sigma_2 = {}$' \
.format(pi1, mu1, s1, mu2, s2), fontsize = 17)
plt.tight_layout()
plt.show()
注意事項與討論：
在樣本數 n = 50 和 n = 100 時，MLE 和 GMM 的估計結果均存在偏差，且兩者之間的趨勢也不一樣。相較之下，GMM 的估計曲線較明顯可看出雙峰形狀，但兩者對於分佈的高峰與尾部的擬合均不理想。
當樣本數增加至 n = 300 和 n = 500 時，MLE 和 GMM 的估計曲線逐步接近真實分佈。此時，MLE 和 GMM 的估計曲線已幾乎重疊在一起。
當樣本數增至 n = 1000 和 n = 5000 時，MLE 和 GMM 的估計曲線與真實分佈曲線十分接近，顯示在大樣本的情況下，兩種方法的估計性能差異趨近於零。
在小樣本條件下，GMM 相較於 MLE，展現出更穩定的估計性能，其分佈曲線能更精準地捕捉混合常態分佈的主要特徵，而 MLE 的估計結果則受樣本波動影響較大，容易導致較大的偏差。
隨著樣本數增加，MLE 的估計性能逐漸提升，與 GMM 的結果趨於一致，兩者在大樣本的情況下均能準確地反映混合常態分佈的結構特徵，顯示出相似的估計能力和趨勢。
本實驗顯示，在混合常態分佈中組成分佈參數的差異較小時，雖然 MLE 和 GMM 在大樣本下均能提供相對準確的估計，但在樣本數不足的情況下，GMM 更具穩定性和準確性，能更有效地捕捉混合常態分佈的主要特徵，因此特別適用於小樣本或資料特徵不明顯的情境。
3.4 評估 MLE 與 GMM 的參數估計表現
與 3.2 的實驗步驟相同，但其中分佈參數改為 \(\pi_1 = 0.3, \mu_1 = -0.5, \sigma_1 = 1, \mu_2 = 1.5, \sigma_2 = 2\)，以此評估 MLE 和 GMM 在參數估計上的表現。
# 設定參數
N = 10000 # number of simulations
n = [50, 100, 300, 500, 1000, 5000]
pi1 = 0.3
mu1, s1 = -0.5, 1
mu2, s2 = 1.5, 2
# 設定演算法初始設定
params0 = [0.3, -0.5, 1, 1.5, 2]
bnd = [(0, 1), (None, None), (None, None), (None, None), (None, None)]
opts = {'disp': False, 'maxiter': 8000}
mean_mini = np.zeros((len(n), len(params0))); mean_gmm = np.zeros((len(n), len(params0)))
bias_mini = np.zeros((len(n), len(params0))); bias_gmm = np.zeros((len(n), len(params0)))
rmse_mini = np.zeros((len(n), len(params0))); rmse_gmm = np.zeros((len(n), len(params0)))
for j in range(len(n)):
# 生成結果矩陣
RESULT_mini = np.zeros((N, 5))
RESULT_gmm = np.zeros((N, 5))
# 生成來自 mix beta distribution 的資料
n1 = binom.rvs(n[j], pi1)
SAMPLE = np.r_[norm.rvs(mu1, s1, size = (n1, N)), \
norm.rvs(mu2, s2, size = (n[j] - n1, N))]
for i in range(N):
X = SAMPLE[:, i]
neg_log_likelihood = lambda x: -np.sum(np.log(x[0] * norm.pdf(X, x[1], x[2]) + \
(1 - x[0]) * norm.pdf(X, x[3], x[4])))
result = minimize(neg_log_likelihood, params0, method = 'Nelder-Mead', \
bounds = bnd, options = opts) # L-BFGS-B
RESULT_mini[i, :] = result.x
# Gaussian Mixture 演算法
gmm = GaussianMixture(n_components = 2, covariance_type = 'full', \
weights_init = [pi1, 1-pi1], \
means_init = np.array([[mu1], [mu2]]).reshape(2, 1), \
verbose = 0, max_iter = 8000, tol = 1e-6)
gmm.fit(X.reshape(-1, 1))
# a, b = np.argmax(gmm.weights_), np.argmin(gmm.weights_)
RESULT_gmm[i, :] = [gmm.weights_[0], gmm.means_[0][0], np.sqrt(gmm.covariances_[0][0][0]), \
gmm.means_[1][0], np.sqrt(gmm.covariances_[1][0][0])]
# minimize
mean_mini[j, :] = np.mean(RESULT_mini, axis = 0)
bias_mini[j, :] = mean_mini[j, :] - params0
rmse_mini[j, :] = np.sqrt(np.mean((RESULT_mini - params0) ** 2, axis = 0))
# gmm
mean_gmm[j, :] = np.mean(RESULT_gmm, axis = 0)
bias_gmm[j, :] = mean_gmm[j, :] - params0
rmse_gmm[j, :] = np.sqrt(np.mean((RESULT_gmm - params0) ** 2, axis = 0))# 整理表格
columns = ['pi1={}'.format(pi1), 'mu1={}'.format(mu1), 's1={}'.format(s1)
, 'mu2={}'.format(mu2), 's2={}'.format(s2)]
# 创建 MultiIndex 列名
columns_mean = pd.MultiIndex.from_product([['Mean for the estimated parameters at N = {}'.format(N)], pd.Index(columns, name='n')])
columns_bias = pd.MultiIndex.from_product([['Bias for the estimated parameters at N = {}'.format(N)], pd.Index(columns, name='n')])
columns_rmse = pd.MultiIndex.from_product([['RMSE for the estimated parameters at N = {}'.format(N)], pd.Index(columns, name='n')])
# 创建新的 MultiIndex
methods = np.repeat(['optimize.minimize', 'GaussianMixture'], len(n))
ns = np.tile(n, 2)
multi_index = pd.MultiIndex.from_arrays([methods, ns])
# 创建 DataFrame 并设置 MultiIndex
df_mean = pd.DataFrame(np.vstack([mean_mini, mean_gmm]), index=multi_index, columns=columns_mean)
df_bias = pd.DataFrame(np.vstack([bias_mini, bias_gmm]), index=multi_index, columns=columns_bias)
df_rmse = pd.DataFrame(np.vstack([rmse_mini, rmse_gmm]), index=multi_index, columns=columns_rmse)
# 使用 Styler 对象来设置列名和单元格内容的对齐方式
df_mean1 = df_mean.style.set_table_styles({
('Mean for the estimated parameters at N = {}'.format(N), 'pi1={}'.format(pi1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Mean for the estimated parameters at N = {}'.format(N), 'mu1={}'.format(mu1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Mean for the estimated parameters at N = {}'.format(N), 's1={}'.format(s1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Mean for the estimated parameters at N = {}'.format(N), 'mu2={}'.format(mu2)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Mean for the estimated parameters at N = {}'.format(N), 's2={}'.format(s2)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
}, overwrite=False).set_properties(**{'text-align': 'center'})
df_bias1 = df_bias.style.set_table_styles({
('Bias for the estimated parameters at N = {}'.format(N), 'pi1={}'.format(pi1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Bias for the estimated parameters at N = {}'.format(N), 'mu1={}'.format(mu1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Bias for the estimated parameters at N = {}'.format(N), 's1={}'.format(s1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Bias for the estimated parameters at N = {}'.format(N), 'mu2={}'.format(mu2)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('Bias for the estimated parameters at N = {}'.format(N), 's2={}'.format(s2)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
}, overwrite=False).set_properties(**{'text-align': 'center'})
df_rmse1 = df_rmse.style.set_table_styles({
('RMSE for the estimated parameters at N = {}'.format(N), 'pi1={}'.format(pi1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('RMSE for the estimated parameters at N = {}'.format(N), 'mu1={}'.format(mu1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('RMSE for the estimated parameters at N = {}'.format(N), 's1={}'.format(s1)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('RMSE for the estimated parameters at N = {}'.format(N), 'mu2={}'.format(mu2)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
('RMSE for the estimated parameters at N = {}'.format(N), 's2={}'.format(s2)): [{'selector': 'th', 'props': [('text-align', 'center')]}],
}, overwrite=False).set_properties(**{'text-align': 'center'})
# 添加水平線在每個 row
df_mean1.set_table_styles([{'selector': 'tr', 'props': [('border-top', '1px solid #666666')]}], overwrite=False)
df_bias1.set_table_styles([{'selector': 'tr', 'props': [('border-top', '1px solid #666666')]}], overwrite=False)
df_rmse1.set_table_styles([{'selector': 'tr', 'props': [('border-top', '1px solid #666666')]}], overwrite=False)
df_mean1.set_table_styles([{'selector': 'tr', 'props': [('border-bottom', '1px solid #666666')]}], overwrite=False)
df_bias1.set_table_styles([{'selector': 'tr', 'props': [('border-bottom', '1px solid #666666')]}], overwrite=False)
df_rmse1.set_table_styles([{'selector': 'tr', 'props': [('border-bottom', '1px solid #666666')]}], overwrite=False)
# 显示 DataFrame
display(df_mean1)
display(df_bias1)
display(df_rmse1)| Mean for the estimated parameters at N = 10000 | ||||||
|---|---|---|---|---|---|---|
| n | pi1=0.3 | mu1=-0.5 | s1=1 | mu2=1.5 | s2=2 | |
| optimize.minimize | 50 | 0.511806 | -0.410114 | 0.932643 | 2.422942 | 1.520930 |
| 100 | 0.457256 | -0.363650 | 0.983922 | 2.289935 | 1.664048 | |
| 300 | 0.398756 | -0.414877 | 1.000286 | 1.872713 | 1.860578 | |
| 500 | 0.338129 | -0.419102 | 1.007015 | 1.782254 | 1.903221 | |
| 1000 | 0.326644 | -0.464029 | 1.008744 | 1.645060 | 1.949919 | |
| 5000 | 0.308993 | -0.497045 | 1.002197 | 1.517628 | 1.994863 | |
| GaussianMixture | 50 | 0.536491 | -0.300080 | 1.023965 | 2.437776 | 1.460377 |
| 100 | 0.484900 | -0.321245 | 1.072371 | 2.320577 | 1.635958 | |
| 300 | 0.403237 | -0.417934 | 1.030275 | 1.855354 | 1.865371 | |
| 500 | 0.345599 | -0.416403 | 1.048865 | 1.778022 | 1.903278 | |
| 1000 | 0.325651 | -0.468830 | 1.022883 | 1.625791 | 1.959839 | |
| 5000 | 0.321356 | -0.484566 | 1.023927 | 1.547544 | 1.991113 | |
| Bias for the estimated parameters at N = 10000 | ||||||
|---|---|---|---|---|---|---|
| n | pi1=0.3 | mu1=-0.5 | s1=1 | mu2=1.5 | s2=2 | |
| optimize.minimize | 50 | 0.211806 | 0.089886 | -0.067357 | 0.922942 | -0.479070 |
| 100 | 0.157256 | 0.136350 | -0.016078 | 0.789935 | -0.335952 | |
| 300 | 0.098756 | 0.085123 | 0.000286 | 0.372713 | -0.139422 | |
| 500 | 0.038129 | 0.080898 | 0.007015 | 0.282254 | -0.096779 | |
| 1000 | 0.026644 | 0.035971 | 0.008744 | 0.145060 | -0.050081 | |
| 5000 | 0.008993 | 0.002955 | 0.002197 | 0.017628 | -0.005137 | |
| GaussianMixture | 50 | 0.236491 | 0.199920 | 0.023965 | 0.937776 | -0.539623 |
| 100 | 0.184900 | 0.178755 | 0.072371 | 0.820577 | -0.364042 | |
| 300 | 0.103237 | 0.082066 | 0.030275 | 0.355354 | -0.134629 | |
| 500 | 0.045599 | 0.083597 | 0.048865 | 0.278022 | -0.096722 | |
| 1000 | 0.025651 | 0.031170 | 0.022883 | 0.125791 | -0.040161 | |
| 5000 | 0.021356 | 0.015434 | 0.023927 | 0.047544 | -0.008887 | |
| RMSE for the estimated parameters at N = 10000 | ||||||
|---|---|---|---|---|---|---|
| n | pi1=0.3 | mu1=-0.5 | s1=1 | mu2=1.5 | s2=2 | |
| optimize.minimize | 50 | 0.333628 | 0.648350 | 0.443891 | 1.640991 | 0.779456 |
| 100 | 0.301684 | 0.615862 | 0.423985 | 1.473734 | 0.613319 | |
| 300 | 0.212067 | 0.329545 | 0.278506 | 0.829515 | 0.318679 | |
| 500 | 0.177408 | 0.313178 | 0.271403 | 0.679992 | 0.245497 | |
| 1000 | 0.121720 | 0.184119 | 0.173141 | 0.426250 | 0.153193 | |
| 5000 | 0.038210 | 0.064305 | 0.065187 | 0.109172 | 0.038963 | |
| GaussianMixture | 50 | 0.333747 | 0.986497 | 0.440079 | 1.584866 | 0.799004 |
| 100 | 0.294395 | 0.704151 | 0.400045 | 1.389733 | 0.617506 | |
| 300 | 0.191937 | 0.304708 | 0.250829 | 0.728294 | 0.307012 | |
| 500 | 0.149357 | 0.271344 | 0.244842 | 0.577402 | 0.235400 | |
| 1000 | 0.094835 | 0.156997 | 0.154378 | 0.328272 | 0.129830 | |
| 5000 | 0.038868 | 0.062473 | 0.065966 | 0.103958 | 0.037234 | |
注意事項與討論：
從結果可以看出，本次實驗的表現與 3.3 有些差異。無論是 MLE 還是 GMM，兩者在平均值、偏差和 RMSE 上的表現相近，難以判定哪種方法更為優越。
此外，兩種方法在樣本數較少時，其估計的平均值與真實值之間存在明顯差異。然而，隨著樣本數增加，平均值逐漸趨於真實值，偏差和 RMSE 也隨之下降。
本次實驗的結果與 3.3 不同，更接近 3.2 的情況。這可能是因為經過 10000 次蒙地卡羅模擬後，兩種方法的估計偏差被進一步縮小，因此難以顯現出 MLE 和 GMM 在參數估計上的顯著差異。
本實驗表明，儘管在實驗次數較少且樣本數不足的情況下，GMM 的表現會優於 MLE，但當樣本數充足或模擬次數足夠多時，兩種方法的參數估計性能幾乎無差異，證實它們都是參數估計中相對穩健的方法。